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Electromagnetic wave scattering from planar dielectric films deposited on one-dimensional, ran- 
domly rough, perfectly conducting substrates is studied by numerical simulations for both p- and 
s-polarization. The reduced Rayleigh equation, which is the integral equation satisfied by the scat- 
tering amplitude after eliminating the fields inside the film, is the starting point for the simulation. 
This equation is solved numerically by considering a random surface of finite length, and by intro- 
ducing wave number cut-offs in the evanescent part of the spectrum. Upon discretization, a system 
of linear equations is obtained, and by solving this matrix system for an ensemble of surface realiza- 
tions, the contribution to the mean differential reflection coefficient from the incoherently scattered 
field, {dRv/d8)- h (v=p,s), is obtained nonperturbatively. It is demonstrated that when the 
scattering geometry supports at least two guided waves, {dRv/d6) iucoh , has, in addition to the well 
known enhanced backscattering peak, well-defined satellite peaks in agreement with theory, for most 
of the parameters considered. 

PACS numbers: 42.25.Fx, 42.25.Gy, 78.66.Bz, 73.20.Mf 
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Over the last five years, the scattering of light from bounded systems where one of the interfaces is rough has 
been studied in several papers. It was demonstrated that these systems may give rise to special enhancements in the 
angular distribution of the intensity of the light scattered incoherently, in addition to the more well-known enhanced 
backscattering peak |l],|^] which is known to appear in the retroreflection direction. These enhancements, known as 
satellite peaks, are present only for bounded systems supporting two (and not many more) surface or guided waves at 
the frequency co of the incident light |J. They are the result of the coherent interference of multiply-scattered waves 
that are time-reversed partners of each other. 

The scattering systems previously considered in the literature in connection with the study of satellite peaks natu- 
rally divide into two classes. The first one consists of a rough dielectric film deposited on a planar perfectly conducting 
substrate [||-|| ■ The second geometry is a thin metal plate whose illuminated (top) surface is a randomly rough one- 
dimensional surface, while the lower one is planar (see also Ref. H). In addition, systems containing volume 
disorder, in contrast to surface disorder, have also been considered, but such systems will not be treated in the present 
work [p. Of . The formalisms applied to the study of these scattering systems range from perturbation theories, such 
as small-amplitude perturbation theory H and many-body perturbation theory MM, stochastic functional meth- 



ods |5|,|9[, and numerical simulations It should be noted that when low-order perturbation theory is not valid, 
a non-perturbative approach has to be taken. Such an approach is provided by numerical simulation techniques like 
the one presented in the present work. 

The scattering geometries considered in the works mentioned above (except Ref. Q) have one thing in common, 
namely that the illuminated surface is always the randomly rough surface, while the back surface is planar. In the 
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present paper we consider a geometry where the lower interface is randomly rough and the upper (illuminated) one 
is planar. In particular, we consider a planar dielectric film deposited on a perfectly conducting rough substrate (see 
Fig. |l|) . Proceeding from the reduced Rayleigh equation for the scattering amplitude we use numerical simulations to 
obtain the contribution to the mean differential reflection coefficient from the incoherent component of the scattered 
light (also known as the mean incoherent differential scattering cross section). 

This paper is organized as follows. We start by describing the scattering geometry, the statistical properties of 
the rough surface, and the satellite peaks supported by this geometry (Sec. ||). In Sec. Ill, the reduced Rayleigh 
equation is derived, and in Sec. |y| the numerical method used in the present work for its solution is discussed. The 
numerical results and their discussion are presented in Sec. [V]. Finally, in Sec. VI the conclusions from these results 
are presented. 



II. THE SCATTERING GEOMETRY 



The scattering geometry we consider in this work consists of a dielectric film, characterized by the dielectric function 
e(uj) = £i{uj) + i£2(w), deposited on a randomly rough perfectly conducting semi-infinite substrate. The top interface 
of the dielectric film is planar, and the medium above the film is assumed for simplicity to be vacuum. The mean 
thickness of the film is d. The scattering system is depicted in Fig. [IJ. The interface separating the perfect conductor 
and the film is described by the surface profile function £(xi). This function is assumed to be a single- valued function 
of x\. Furthermore, it will be assumed to be a zero- mean, stationary, Gaussian random process defined through the 
following properties 

(C(X!))=0, (2.1) 
(C(zi)C(z'i)) =<J 2 W{\x 1 -x' 1 \). (2.2) 

Here a is the RMS-height of the surface roughness, and is the surface-height autocorrelation function nor- 

malized such that W(0) = 1. The angle brackets denote an average over the ensemble of realizations of the surface 
profile function C( x i)- It will later on prove useful to have the power spectrum, g(|fc|), of the surface roughness at our 
disposal. It is defined as 

/oo 
dx 1 W{\x l \)e- ikx K (2.3) 
-oo 

In the present paper two forms of the power spectrum will be considered. These are the Gaussian power spectrum 
defined by 

g(\k\) = V^ae-^ ) (2.4) 
where a is the transverse correlation length of the surface roughness, and the West-O'Donnell power spectrum 

g (\k\) = - [0(k - k_)6{k + -k) + 0(-fc_ - k)d{k + k+)] , (2.5) 
k + — fc_ 

where 6{k) is the Heaviside unit step function. The latter power spectrum was recently used by West and O'Donnell |l2] ] 
in an experimental study of enhanced backscattering by the surface plasmon polariton mechanism in the scattering of 
p-polarized light from a one-dimensional, randomly rough, metal surface. In their work West and O'Donnell defined 
k + and fc_ by 

k± = k sp (uj) ± {lo/c) sin6> wo , (2.6) 

where ±fc sp (w) are the wave numbers of the forward- and backward-propagating surface plasmon polaritons at a 
planar vacuum-metal interface, whose frequency to is that of the incident field. The physical significance of the angle 
9 wo is that if the angle of incidence Oq is in the interval (— 9 wo , WO ), the incident light can excite both the forward and 
backward propagating surface plasmon polaritons through the surface roughness. Similarly, if the scattering angle 
9 is in the interval (—0 WO ,0 WO ), the excited surface plasmon polaritons will be coupled to scattered volume waves 
in the vacuum region above the metal surface. With this form of the power spectrum the incident light couples 
strongly to the surface plasmon polaritons over a limited range of angles of incidence rather than weakly over a large 
range of this angle, as is the case when a power spectrum g(|fc|) peaked at k = 0, e.g. a Gaussian, is used. This 
is because in the latter case the wave numbers ±k sp (uj) lie in the wings of the power spectrum, where it is usually 
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small. As a co nsequence of the strong excitation of surface plasmon polaritons by the incident light when the power 
spectrum (|2.5| ) is used, the amplitude of the enhanced peak, which is caused by the coherent interference of multiply 
scattered surface plasmon polaritons with the ir time-reversed partners, is significantly increased with respect to its 



value when the Gaussian power spectrum (2.4) is used. 

To get an intuitive picture of these surfaces it is useful to supply the mean slope, s, and the mean distance between 
consecutive peaks and valleys, (D), as measured along the (lateral) a;i-direction. For a stationary zero- mean, Gaussian 
random process, the RMS-slope s is related to the power spectrum by 



and a good estimator for (D) has been shown to be [H 



(D) 



k*g(\k\)' 

For the two power spectra considered in this work, these quantities are given by 

{ y/2— Gaussian 



and 



-^\Jk 2 + + k+k- + k 2 _ West-O'Donnell 



-?=a Gaussian 

V6 



(2.9) 



L>! ' ^/tS—IJ West-O'Donnell ' 



Two surface profiles with the same (Gaussian) height distribution, but with a Gaussian and a West-O'Donnell power 
spectrum possessing nearly the same value of the RMS-slope s are plotted in Fig. 0. 

A. Satellite peaks 

Satellite peaks appearing in the scattering of light from bounded systems are well-defined enhancements in the 
angular distribution of the scattered (or transmitted) intensity evenly distributed around the enhanced backscattering 
peak. They are multiple-scattering phenomena, and are a consequence of the constructive interference of multiply- 
scattered waves with their time-reversed partners. In order for satellite peaks to exist, the scattering system must 
support N guided waves, where N is at least two. If we denote the wavenumbers of these N guided waves by 
qi(uS) ,q2(ui) , . . ., qpf(w), it can be demonstrated that in the absence of roughness and absorption (es(tjj) — 0), the 
satellite peaks should appear at scattering angles defined by the following relation 

s[n0 t"i,n) = - s ' m ®o ± - [<7mM - gn (<*>)] , m^n. (2.11) 
It should be stressed that not all of these angles may correspond to real satellite peaks. Some of them may lie in the 



evanescent (non-radiative) part of the spectrum, i.e. the right hand side of Eq. (|2.11) may be greater than unity in 



magnitude. Furthermore, not all of the real satellite peaks are guaranteed to be observed, due to their low intensity. 

In the present work when the power spectrum ( |2.5| ) is assumed, we will choose the values of k + and fc_ to include 
the wavenumbers of selected guided waves supported by the scattering system depicted in Fig. 1, enhancing the 
roughness-induced excitation of these guided waves as a result. We will see that the amplitudes of the corresponding 
satellite peaks, as well as the amplitude of the enhanced backscattering peak, will be increased thereby, in comparison 
with their values when a Gaussian power spectrum is used. 

It should be obvious that the wavenumbers of the guided waves, q n (io), are related to the mean thickness, d, of the 
dielectric film. By applying the wave equation with the proper boundary conditions, the dispersion relation for the 
guided waves can be obtained. In the case of no roughness, i.e. C( x i) = 0, and no absorption, it is given by 

e((j)0o(q,(j) = a(q,cj) tan a(q, to)d (2-12) 
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for p-polarization p|, and by 

/?o(g, lo) — a(q, u>) cot a(q, to)d (2-13) 
for s-polarization 0], with (3o(q,uj) = y/ q 2 — (w 2 /c 2 ) and 



a(q, uj) — de(uj)^j — q 2 , Rea(q,cu) > 0, Ima(q,uj) > 0. (2-14) 

From these dispersion relations, the number of guided waves supported by the scattering geometry can be deduced. 
This number n is controlled by the following inequalities for p-polarization ^ 

n — 1 d n „. 

— ; < t < — ; , 2.15) 

2y/F^l A 2^/7^1 

and for s-polarization j| 

2n - 1 d 2n + 1 . 

; < — < ; 2.16 

where n = 1,2,3,... is the number of guided modes supported by the system at the wavelength A of the incident 
light. 

III. THE REDUCED RAYLEIGH EQUATION 

In the present work we will consider a plane wave incident from the vacuum side on the structure described 
earlier (Fig. |l|), with the xia^-plane the plane of incidence. By this choice for the plane of incidence, we are guaranteed 
that there is no cross-polarized scattering, and the plane of scattering coincides with the plane of incidence. As a 
consequence, one field component is sufficient to describe the electromagnetic field completely. Thus, to simplify the 
notation we introduce 

= $ v (x u x 3 \Lj)e- iut e 2 , (3.2) 

where e 2 is a unit vector in the a^-direction, and the subscript v is a polarization index (y = p, s). We have explicitly 
assumed a time-harmonic dependence of the fields here. 

In the region above the film, 2:3 > d, the field consists of an incident and scattered waves. It can be represented by 

<f> u (x 1 ,x 3 \co) = e ikxl - iao( - k ' oj)x3 + f ^R v (q\k)e iqXl+iao ^ X3 , (3.3) 



where R v (q\k) is the scattering amplitude, and 



q 2 , \q\ < 7. 



a (q,u)=i V? 1_ - (3.4) 



^q 2 -^, \Q\>*. 

Furthermore, the field inside the film consists of both upward and downward propagating modes. By applying the 
boundary conditions satisfied by the field at the vacuum-dielectric and the dielectric-metal interfaces, two coupled 
integral equations are obtained (the Rayleigh equations). In obtaining these equations the Rayleigh hypothesis has 
been imposed 0,^5). The condition for the validity of this hypothesis can crudely be stated as ]l6) <C 1. 

Now, by eliminating the modes inside the film, a single integral equation for R v {q\k), known as the reduced Rayleigh 
equation, is obtained. For the scattering system considered in the present work it assumes the following form |17J 

^M+(p\q)R v (q\k) = M~{p\k), (3.5) 

where 



4 



PI 



-ia( g ,u,)d a(q,u)±e(uj)ao{q,uj) 



I(a(q,cj)\p- q) 



+ eiaMd «Mi^ ; ( _ a( w)| _ j 



for p-polarization, and 



M±(?|p)=±e 



<(4iW)d /(-gfoaQjp-g) 



for s-polarization. The function I(j\q) appearing in these formulae is defined by 

poo 



(3.6) 



(3.7) 



(3.8) 



and a(q,cu) was defined earlier in Eq. (2.14). The main purpose of this paper is to calculate the scattering amplitude 
R v (q\k) numerically. This amplitude is related to a physically measurable quantity, the mean differential reflection 
coefficient (dR v /d6), which is defined as the fraction of the total incident flux that is scattered into a small angular 
interval d6 around the scattering direction 9. Since this quantity for weakly rough surfaces will have a dominating 
peak due to coherent (specular) reflection, it is useful to subtract off this contribution. If we do so, we are left with 
what is called the contribution to the mean differential reflection coefficient from the incoherent component of the 
scattered light, which we will denote by (dR u /dO), ncoh . This quantity is related to the scattering amplitude R v (q\k) 
according to M 



8R V 
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incoh 



1 



L 27rc cos#o 



(\R„n-\(R u ) 



(3.9) 



Here L is the length of the sample along the xi-direction, and #0 and 9 are the angles of incidence and scattering, 
respectively. These angles are defined in the counter and clockwise directions as indicated in Fig. [I], and they are 
related to the wave numbers k and q by 



sin t>a, 



sin0. 



(3.10) 



IV. THE NUMERICAL METHOD 

In the paper by Madrazo and Maradudin [O , their reduced Rayleigh equation was solved numerically by replicating 
the rough surface of length L an infinite number of times. By doing so, they covered the entire a; 1 -ax is and obtained 
a diffraction grating of period L. Hence the wave number integration in their equivalent of Eq. ([T^) was converted 
into an infinite, but for practical implementation a large but finite, sum over Bragg beams. By this method they were 
able to obtain convergent results for (dR v j '89) i nco h- 

In the present paper, we have chosen probably the most straightforward method possible for the numerical solution 
of the reduced Rayl eigh equation. We do not replicate the surface periodically, but instead truncate the wave number 



integration in Eq. (3^5) somewhere in the evanescent (non-radiative) part of the spectrum, say at q — ±Q/2, where 
Q ^> uj/c. With a finite length of the surface L ^> A, the reduced Rayleigh equation can now be discretized by a 
standard quadrature scheme, and a system of linear equations is obtained. The result is 

h Nq ' 2 N N 

^ V w n M+(p m \q n )R v ( qn \k) = M-(p m \k), (4.1) 

n=-N q /2 

where N q + 1 is the number of discretization points in wave number space, h q — Q/N q is the corresponding discretiza- 
tion length, and {w n } are the weights of the quadrature scheme used. Furthermore, the abscissas {q n } arc defined 
by 
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q n = nh q , n = -^,...,^-. (4.2) 

The quadrature scheme used in the numerical results presented in the next section is an 0(h*) method with @ 
wi = u>jv = 3/8, W2 = wn- i = 7/6, = wn-2 = 23 /24, and with all other weights equal to one. The matrix 
elements M^(p m \q n ) in Eq. (|4.l|) , are given by Eqs. (|3.6| ) and (3/7), but now with a finite length of the surface, viz. 



J( 7 | ? )= / ^dxie l ^ {xi) e- tqxi . (4.3) 
Integrals of this f orm ar e o ften referred to as Fourier integrals. With finite wave number cut-offs at q — ±Q/2, it 



follows from Eqs. (3.5)— (3.7) that wave numbers in the range [-Q, Q] are needed in the calculation of these integrals. 
Since, as we will see in the following subsection, these integrals will be evaluated by discretizing the spatial a^-variable 
and taking advantage of the (discrete) Fourier transform, the value of Q is controlled by the spatial discretization 
length h of the problem. To resolve wave numbers up to Q, one has from the relation for critical sampling (Nyquist 
freq uenc y) that the number of spatial discretization points has to satisfy n/h > Q = h q N q , where we have used 



Eq. (|4.2j ). In principle the wave number discretization length used in Eq. (4.1), h q , and the one obtained from the 
Fourier transform, h' q , are independent. However, from a numerical point of view this is not very practical, and we 
will chose them to be the same, h' q = h q . From the theory of the Fourier transform we know that the discretization 
lengths in real and wave number space are related by h q — 2ir/(Nh), where N is the number of spatial discretization 
points. If this latter expression is used in n/h > Q = h q N q , one arrives at 

N > 2N q . (4.4) 

In the numerical calculation we have chosen N = 2N q in order to avoid unnecessary calculations. 



A. The integrals I(y\q) 



In order to calculate the Fourier integrals in Eq. (4.3), care has to be taken. The reason for this is the oscillatory 
nature of these integrals which appears unavoidable for large absolute values of the wave number parameter q. For 
example, if one use a standard discretization approach to this integral, I{^\q n ) = h ■ exp(iq n L/2) Y^, m =o ^\ e ^P(^Cm)], 
where T denotes the Fourier transform, one will likely obtain inaccurate and even wrong results (see Ref. [Q Sec. 13.9 
for details). In a more sophisticated and reliable method, the integrand is approximated by higher-order (piecewise La- 
grange) interpolating polynomials. By taking advantage of the additional information provided by these interpolating 
polynomials, highly accurate results for the oscillating Fourier-integrals can be obtained from 



I(j\q n ) ~ h 



W(q n h)T[{e^}] 



(4.5) 



whe re th e j'-summation extends over the endpoints of the integration interval. This endpoint correction appearing in 
Eq. (4.5) enters due to the difference in the functional form of the interpolating polynomials for internal and boundary 
points. The specific form of W(-) and Tji depends only on the order of the interpolating polynomials. In this work 
cubic order has been used, and the interested reader is referred to Ref. Jl8| for explicit expressions for these functions. 

The most time consuming part of the numerical solution of the reduced Rayleigh equation is not, as one might 
expect at first glance, the solution of the matrix system. Instead it is to set up this system and, in particular, to 
calculate the integrals I{p(\q). We stress that these integrals in principal have to be calculated for each new value of q. 
This means, for the method just sketched, that the Fourier transform has to be recalculated for each new value of the 
wave number variable q. Even though the Fourier transform has an fast implementation, the fast Fourier transform 
(FFT), there are so many transforms to be calculated that the whole operation becomes quite time consuming. For 
a system of linear size N, the number of integral evaluations needed is 2N 2 , where the factor of 2 enters because one 
needs two integral evaluations per matrix element. Another method, which overcomes this drawback when applic able , 
was introduced in Ref. |fTT|| , and consists of expanding the exponential function, exp(i^Q{xi)) appearing in Eq. (4.3) 
in a Taylor series, and integrating the resulting expression term-by-tcrm 



I(7\Q) 



E 

n=0 



dx\ 



(»7) r 



OO 

E 

71 = 



(n) r 



C n (zi) 



(4.6) 
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For surface profiles of modest roughness, this expansi on w ill converge rather quickly, meaning that only a limited 



number of terms are needed in the summation in Eq. (4J3). However, the real advantage of this method is that the 



Fourier-transforms needed to calculate I{^\q) do not have to be recalculated for every new value of the wavenumber 
variable 7, in contrast to what is the case for the Fourier-method presented above. All this, as we will see, results in 
a dramatic reduction of computational time (cf. Fig. ^|). 

For the roughness used in the present work, the latter method for the calculation of I("f\q) will be the method 
of choice. The Fourier-method has been included and used to demonstrate and compare the results obtained by 
its use to those obtained by the Taylor-method. However, the Fourier-method is the method of choice when the 
Taylor-method begins to converge slowly, which can occur for large-amplitude roughness. When we are forced to use 
the Fourier-method, high demands are placed on computer time. 

V. RESULTS AND DISCUSSIONS 

In this section we present the results of the simulations obtained by the method just described. For all simulations 
incident light of wavelength A = 633 nm was used. The rough surface had length L — 160. 1A = 101343.3 nm, and 
an RSM-height a — 30 nm, unless indicated otherwise. The rough surfaces were generated by the method described 
in Refs. |^|n|- The number of surface realizations used in the ensemble average was = 3000, if nothing is said 
to the contrary. Such a large number of samples was used in order to reduce the noise level which such numerical 
calculations are plagued with. If nothing else is said, the integration method used in the calculation of the integrals 
/(7I?) was the Taylor-method where ten terms were retained. This was enough, as we will see, to obtain convergent 
results. The calculations about to be presented were performed on an SGI/Cray Orion 2000 supercomputer, and the 
typical CPU-time for this number of samples was roughly 30 hours. Furthermore, the number of spatial discretization 
points was set to TV = 1604 (and hence N q = 802). Since only a limited number of Fourier transforms was needed per 
sample when using the Taylor- met hod, the advantage of the FFT-algorithm is rather marginal as compared to the 
system size. 

In the absence of absorption, i.e. when the dielectric function is real, the following unitarity condition (conservation 
of energy) , coming from the unitarity of the scattering matrix, should be satisfied 

p- " o(g,h,) R*(q\k)R(q\p)=2 7 rS(k-p). (5.1) 

Numerical simulations with negligible absorption showed that this condition was satisfied within 0.5% in the case that 
k = p. With a small non- vanishing imaginary part of the dielectric function (£2(^1) = 0.01), this relation was satisfied 
within 4%-16%, depending on the roughness used. 

A. S-polarization 

We start the presentation of the numerical results by considering s-polarized incident light. The dielectric constant 
of the film is, at the given wavelength, e(uj) = 2.6896 -M0. 01, and the mean film-thickness is set to d = 500 nm. This is 
the same set of parameters used for s-pol arizat ion in Ref. pT|. In the absence of roughness and absorption, this mean 



film thickness predicts, according to Eq. (2.13), two guided wave modes with wave numbers qi(u>) — 1.5466(w/c) and 
<7 2 (w) = 1.12423(cj/c) (cf. Ref. Qj). The corresponding satellite peaks are then, according to Eq. ( p. 11 ), expected to 
appear at 9 = ±17.7°. 

In Fig. ^| we present the results of the numerical simulation of the scattering of s-polarized light normally incident 
on a surface for which both the surface height distribution and the surface-height autocorrelation function have the 
Gaussian form. The correlation length of the surface roughness is a = 100 nm. With these parameters the RMS-slope 
and mean distance between consecutive peaks and valleys are given by s = 0.424 and (D) = 128.3 nm, respectively. 
In the raw data for the mean incoherent reflection coefficient (Fig. |^a) the satellite peaks are hard to see due to the 
noisy background, but the enhanced backscattering peak at 6 — 0° is easily located. This conclusion is not affected 
by doubling the number of Taylor-terms retained in the calculation (result not shown). However, we believe that 
the noise stems from the fact that we are using a plane incident wave instead of a beam of finite width. In order to 
reduce the effect of the noise, we have applied a standard five-point filter to the raw data. This filter affects only 
eleven consecutive points, and it should be noted that this smoothing procedure is not as aggressive as the one used 
by Madrazo and Maradudin [iljj . After applying the smoothing-filter to the raw data, the satellite peaks appear 
more clearly (Fig. ^b). Even though their amplitudes are quite small, they seem to appear at the correct angles as 
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indicated by the dashed line in the figure. Satellite peaks, like enhanced backscattering peaks, are multiple scattering 
effects, and can consequently be masked by single-scattering phenomena. This is probably the reason for the small 
amplitudes of the satellite peaks in Fig. ||b. 

We now focus on the West-O'Donnell power spectrum. From perturbation theory ||, it is known that the lowest 
order contribution, i.e. the single-scattering contribution, to the mean incoherent differential reflection coefficient is 
proportional to the power sp ectru m g(\q — k\), where q and k are the wave numbers of the scattered and incident 
waves, respectively (cf. Eq. ( |3.10 )). This implies that the single-scattering contribution to {dR v /d6)- lncoh is largest 
for angles where g(\q — k\) is large. As mentioned earlier, West and O'Donnell constructed a power spectrum g(|fc|) 
which is non- vanishing only in a limited range of |fc| defined by the upper and lower limits k+ and k^, respectively. 
Thus, for angles of incidence and scattering satisfying \q — k\ < /c_, single scattering processes will not contribute 
to {dR„/ ^96*)incoh: an d may enhance the roughness-induced coupling of electromagnetic waves to surface plasmon 
polaritons if the values of k± are properly chosen. In the scattering process we are about to consider, we have chosen 
fc_ = 0.82(w/c) and fc+ = 1.97(w/c), and thus s = 0.427 and (D) = 201.1 nm, so that the incident wave, for small 
angles of incidence, can excite both guided waves supported by the geometry through the surface roughness. With 
these parameters, and for normal incidence, one finds from Eq. ( 3.10 ) that single-scattering effects do not contribute 
to the mean incoherent differential reflection coefficient for scattering angles \8\ < 55.1°. Since the satellite peaks, 
located at 9 = ±17.7° in this case, cannot be masked by single-scattering processes, one now expects more pronounced 
satellite peaks then those obtained with the use of the Gaussian correlation function presented in Fig. 3. As can be 
seen from Fig. |Ja, this is indeed what we get. The satellite peaks are well separated from the background in the raw 
data, and no smoothing is needed. The amplitude of the prominent enhanced backscattering peak seen in this fi gur e 
is twice that of its background. This is another indication that single scattering processes have not contributed [g0| . 
The abrupt increase in (dR s /d9) incoh for angles \9\ somewhat above 50° is due to single-scattering effects setting in. 
We predicted above that this should happen for angles \9\ > 55.1°, a result that fits quite well with what can be 
read off from the figure. In Fig. [|b we have shown the (smoothed) result of the simulation obtained by Madrazo and 
Maradudin jll[] for the same set of parameters, but with the top interface being the rough one instead of the lower 
one. It is observed that the overall amplitude is three times smaller than the one presented in Fig. ^a. Furthermore, 
the intensity obtained for our geometry in the single scattering regime (\9\ > 55.1°) has a broader distribution than 
the one obtained by Madrazo and Maradudin Jll[]. The reason for this has to do with the increased reflection at the 
rough surface present in our case. 

In the previous section we discussed two methods for the calculation of the integrals /(7I5). We will now compare 
these two methods with respect to numerical performance and accuracy. We will use the West-O'Donnell power 
spectrum with the parameters used above, but with the number of spatial discretization points chosen to be N = 
2 11 = 2048. This is done in order to take advantage of the fast Fourier transform. In Fig. || the results for the 
mean incoherent reflection coefficients, (dR s /d9) incoh are presented for the Fourier- (Fig. ||a) and Taylor- (Fig. ||b) 
methods. The number of samples used in obt aining these results was N{ = 50. For the Taylor-method, the number 
of terms retained in the Taylor expansion (4.6) was ten. The CPU-time needed to obtain these results on a SGI/Cray 
Orion 2000 supercomputer were ty = 63.3 CPU- hours and tx = 1.3 CPU-hours for the Fourier- and Taylor-methods, 
respectively. In Fig. ^c the discrepancy between the results of the two methods are shown. The relative error is 
roughly 1%, which is of the order of the error due to the use of a finite number of samples [jl9|| . Thus we may conclude 
that for the parameters used, the results are equivalent, and a factor of 50 in performance is gained by using the 
Taylor-method instead of the Fourier-method. Consequently, the Taylor-method will be the method of choice for 
slightly rough surfaces, and we will use this method in the following numerical calculations. However, we should 
stress that this result is not a general one, but depends heavily on the roughness parameters used. In particular, for 
sufficiently large RMS-heights, the Taylor-method will converge slowly, and we are left with the numerically demanding 
Fourier- method. 

As the angle of incidence departs from the direction of normal incidence (60 = 0°), the amplitudes of the satellite 
peaks are known to decrease. In Fig. ^ we show the results of a simulati on for a geometry with the same parameters 
as above, but now with angle of incidence of 9o — 5°. According to Eq. ( 2.11 ) the satellite peaks should now appear 
at 9 = 12.6° and 9 = —26.1°. It is seen from the figure that the amplitudes of these peaks, and the enhanced 
backscattering peak (at 9 = 5°), are reduced in amplitude compared to the case of normal incidence, but they are 
still easily distinguished from the background and are located at the predicted positions. 



B. p-polarization 



We will now focus our attention on p-polarization of the incident light. Satellite peaks are typically harder to 
observe in p- than in s-polarization. This is related to the reduced reflectivity at the rough surface of p-polarized 
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light as compared to that of s-polarization for the same angle of incidence. In fact, this was the reason why Madrazo 
and Maradudin |p"l|| used different parameters for the two polarizations considered. However, in this paper, the 
rough surface is perfectly reflecting, and we thus use the same set of parameters for the two polarizations. For 
completeness and comparison, at the end of this section we have included results for the parameters used by Madrazo 
and Maradudin Q . 

W e star t by presenting the numerical results obtained for the parameters used in the preceding subsection. From 
Eq. ( [2.15 ) one finds that the scatterin g sys tem, now in the case of p-polarization, supports three guided waves, which 
according to the dispersion relation ( 2.12 ) have the wave-numbers qi(oS) = 1.6125(w/c), q2(u) = 1.3821(w/c) and 
q3(uj) — 1.0029(w/c). He nce there are six possible satellite peaks, and for normal incidence they are located at 
scattering angles (cf. Eq. ([Hi])) 0± 2) = ±13.3°, 0f 2 3) = ±22.3°, and 6f 13) = ±37.6°. 

In Fig. ^ the contribution (dR p /d9)i nco h is presented for a rough surface described by a Gaussian power spectrum 
with correlation length a = 100 ran. As for the case of s-polarization with the same set of parameters, the enhanced 
backscattering peak is easily seen in the raw data (Fig 0a). However, there is now clear evidence of the satellite peaks 
in these data. After smoothing (Fig. ^|b), the satellite peaks at 6% 3 ^ are clearly visible, while there is only a weak 

indication of those at 0^ 2 y However, there is no indication of the satellite peaks at 0^[ neither in the raw, nor in 
the smoothed data. 

We saw earlier that, in going to a surface roughness described by the West-O'Donnell power spectrum, the satellite 
peaks became much more pronounced because they were not masked by single-scattering processes. Figure || presents 
the results for the scattering of normally (9q = 0°) incident p-polarized light when the power spectrum is of the West- 
O'Donnell type with parameters fc_ = 0.82(cj/c) and k + = 1.97(cj/c) (as above). Here, as in the case of s-polarization, 
the satellite peaks at 0^j ^ are well defined both in the raw (Fig. ||a) and smoothed data (Fig. IJd). In the raw data, the 

peaks at 6*^ ^ can also be seen, even though they appear more clearly in the smoothed data. As in the Gaussian case, 
there is no sign of the satellite peaks at 9^ 3 j neither in the raw nor in the smoothed data. This may indicate that the 
coupling of the incident wave to the guided mode with wavenumber 93(0-") is quite weak. We observe that the dominant 
satellite peaks (0^ ^ ) have roughly the same amplitude as those observed for s-polarization. This has to do with the 
fact, as we claimed above, that the randomly rough surface is perfectly reflecting, and thus no distinction between 
p- or s-polarization when it comes to reflectivity has to be made. However, for the geometry considered by Madrazo 
and Maradudin jllj , where the rough surface was the vacuum-dielectric interface, the amplitude of the satellite peaks 
for p-polarization were much lower then those obtained in s-polarization, due to the decreased reflectivity for the 
p-polarized light as compared to that for s-polarized light. 



If the incident light has an angle of incidence of, say, 9q = 5°, the satellite peaks should appear at 0q 



(1,2) 



8.2°, -18.5°, 9% 3) = 17.0°, -27.8°, and 0^ 3) = 31.5°, -44.2°, according to Eq. ( p.ll| ). The results for this situation, 
still with the West-O'Donnell power spectrum, are given in Fig. ^a and|^b for the raw and smoothed data, respectively. 
Also in this case the peaks at 0^ ^ and 0^ ^ are detectable. However, it is interesting to note that the peaks at 0^ 3 %, 
of which there was no sign for an angle of incidence 9q = 0° , are weakly visible in both the raw and smoothed data. 
This we find rather surprising, since typically the satellite peak amplitudes decrease quite rapidly with increasing 
angle of incidence. 

In the Madrazo and Maradudin paper they used another dielectric function, e(uj), and mean film thickness, 
d, for p-polarization than they (and we) used in the simulations for s-polarization. In particular, they used e(lu) = 
5.6644 + O.Oli and d — 380 nm. All other parameters were the ones used earlier in this paper. For the purpose 
of comparison, we have included simulation results for these parameters. In Fig. 10 the result for (dR p /d9) in h 
is presented for a surface possessing a Gaussian power spectrum with a correlation length a — 100 nm (and with 
= 30 nm as before). This corresponds to an RMS-slope of s = 0.424, and the mean distance between consecutive 
peaks and valleys is (D) = 128.3 nm. For such a surface, in the limit of no roughness and no absorption in the film 
(£2(w) = 0), the structure supports three guided modes whose wave numbers are q \{oS) = 2.34(w/c) (72(1-^) = 2.04(w/c) 
and 93 (w) = 1.32(lj/c). Hence, there are six satellite peaks, which according to Eq. ( 2.11 ), are located at 0(i,2) — ±17.5, 
$(2,3) = ±46.1, while 0(i,3) lie in the evanescent (non-radiative) part of the spectrum and thus are not visible. From 
Fig. [To] it is seen that even the satellite peaks in the radiative part of the spectrum are not possible to locate in the 
numerical results. This is in fact the same conclusion drawn by Madrazo and Maradudin JTT[ ] for their geometry. 
Thus the increased reflectivity we have at the rough surface does not seem to affect the detection of the possible 
satellite peaks. In Fig. |TT| results for the West-O'Donnell power spectrum defined by the parameters fc_ = 1.61(a->/c) 
and k+ = 2.76(w/c) (which gives s = 0.658 and (D) = 137.3 nm) are presented. In this case the satellite peaks are 
predicted to appear at angles ± = ±17.5° for normal incidence (0o = 0°) and at ± = 12.3°, —22.8° for 9g = 5°. It 
should be noted that we only expect one pair of satellite peaks because the incident wave can directly excite only two 
of the three possible guided waves (fc_ > 93(0;)). Furthermore, with , k± > uj/c, single-scattering processes will not 







contribute at all (see the discussion above). From Fig. [0] we see that the satellite peaks in the smoothed data are 
found at the expected positions both for 6q — 0° and 9q = 5°, even though their amplitudes are small. This was the 
the same conclusion arrived at in Ref. [O. 



VI. CONCLUSIONS 



We have presented a numerical study of the scattering of electromagnetic waves of both p- and s-polarization 
from a system consisting of a planar dielectric film deposited on a randomly rough perfectly conducting substrate. 
The numerical calculation was performed by considering a finite length of the randomly rough surface L> A, then 
by solving the corresponding reduced Rayleigh equation for the scattering geometry by standard techniques. By 
averaging the results for the scattering amplitude and its squared modulus over the ensemble of realizations of the 
rough surface, the mean incoherent reflection coefficient was obtained. The numerical results for the mean differential 
reflection coefficient shows that the scattering geometry under investigation gives rise to satellite peaks at well defined 
positions for most of the scattering and roughness parameters considered, in accordance with theory 

From a methodological standpoint the present work demonstrates that a purely numerical approach to the solution 
of a reduced Rayleigh equation by standard numerical techniques, without the necessity of replicating a segment of 
random surface of length L periodically, as was done in pd| ] , is a viable, nonperturbative approach to the investigation 
of the scattering of light from one-dimensional random surfaces whose roughness is such that the Rayleigh hypothesis 
is valid. 

From a physical standpoint the results presented in Fig. [| strongly suggest that the scattering system studied in 
the present work is a more favorable one for the experimental observation of satellite peaks than the one studied 
in |^,^ 11 1, due to the increased incoherent scattering intensity to which it gives rise. In addition, structures of the 
kind studied in the present work should be easier to fabricate than the one considered in ^,^,|ll[. 
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FIG. 1. The scattering geometry considered in the present paper. The medium above the illuminated planar surface at x% = d 
is vacuum. The randomly rough surface at x% = C,{x\) is characterized by a stationary, zero-mean, Gaussian random process 
Q(xi). Below this surface the medium is a perfect conductor, while a dielectric film is present in the region (,(x\) < X3 < d. 
This film is characterized by the frequency dependent dielectric function e(ui). The angles of incidence and scattering are 9q 
and 9, as indicated in the figure. 

FIG. 2. Two rough profiles both with a Gaussian height distribution with an RMS-value a = 30 nm. The power spectrum is 
of the (a) Gaussian type, with a — 100 nm, and the (b) West-O'Donnell type, with fc_ = 0.82(oj/c) and k + — 1.97(w/c). Here 
the wavelength is A = 633 nm. With these parameters the RMS-slope and distance between consecutive peaks and valleys are 
respectively s = 0.424, (D) = 128.3 nm, and s — 0.427, (D) = 201.1 nm for the case of the Gaussian and West-O'Donell power 
spectra. Note that there are different scales on the first and second axes, with the result that the profiles appear much rougher 
than they really are. 

FIG. 3. The contribution to the mean differential reflection coefficient from the incoherent component of the scattered light 
{dR s /d8) incoh as a function of the scattering angle 9 when an s-polarized plane wave of wavelength A = 633 nm is incident 
normally (9o = 0°) on the scattering system depicted in Fig. [IJ. The dielectric function of the film is e(lo) = 2.6896+iO.Ol, and its 
mean thickness is d — 500 nm. The surface profile function £(xi) is characterized by a Gaussian surface height distribution with 
an RMS-height, a = 30 nm, and a Gaussian surface-height autocorrelation function with a 1/e-correlation length of a — 100 nm. 
The scattering system supports satellite peaks at ±17.7° as indicated by the vertical dashed lines in the figure. The raw data 
are presented in Fig. |^a, while the their (five-point filter) smoothed analogs are given in Fig. []b. 

FIG. 4. The contribution to the mean differential reflection coefficient from the incoherent component of the scattered light 
{9R s /d9) lncoh as a function of the scattering angle 9 when an s-polarized plane wave of wavelength A = 633 nm is incident 
normally (Bo — 0°) on the scattering system depicted in Fig. The dielectric function of the film is s(u>) = 2.6896 + i0.01, and 
its mean thickness is d — 500 nm. The surface roughness profile function is characterized by a Gaussian surface height 

distribution with an RMS-height a — 30 nm and a West-O'Donnell power spectrum defined by the parameters fc_ = 0.82(oj/c) 
and k+ = 1.97(a>/c). The scattering system supports satellite peaks at ±17.7° as indicated by the vertical dashed lines in the 
figure. The raw data for for the geometry considered are presented in Fig. |j|a. The numerical results obtained by Madrazo and 
Maradudin jll| for a corresponding geometry with the vacuum-dielectric interface being the rough one (see text) are given in 
Fig. |b. 

FIG. 5. The same as Fig. ^a (raw data), but now using the (a) Fourier- and the (b) Taylor-methods for calculating I("y\q) 
(see text). In Fig. the difference between the results for the contributions to the mean differential reflection coefficient 
calculated by these two methods, A (9R s /d9) lncoh , is presented. We observe that the relative discrepancy is roughly of the 
order of 1% for the parameters used. The number of samples used in obtaining the results was N( = 50, and the number of 
spatial discretization points was N = 2 11 = 2048 in order to take advantage of the Fast Fourier transform. For the two methods 
the same surface profiles have been used in the calculation. The other parameters are the one used in obtaining the results 
plotted in Fig. ^. 

FIG. 6. The same as Fig. (raw data), but now for an angle of incidence 6o = 5°. The satellite peaks should now be present 
at scattering angles 9 = 12.1° and 9 = —26.1°, as indicated by the vertical dashed lines. 

FIG. 7. The contribution to the mean differential reflection coefficient from the incoherent component of the scattered light, 
{dR P /d9) incoh as a function of the scattering angle 9 when a p-polarized plane wave of wavelength A = 633 nm is incident 
normally (9o — 0°) on the scattering system depicted in Fig. [j]. The dielectric function of the film is e(ui) = 2.6896 + iO.Ol, 
and its mean thickness is d = 500 nm. The surface profile function C,(x\), is characterized by a Gaussian surface height 
distribution with an RMS- height, a — 30 nm and a Gaussian surface- height autocorrelation function with a 1/e-correlation 
length of a — 100 nm. The scattering system may give raise to six satellite peaks at #^2) = ±13.3°, 0q,3) ~ ±22.3° and 
9^ 3 j = ±37.6° as indicated by the vertical dashed lines in the figure. The raw data are presented in Fig. ^ja, while the their 
(five-point filtered) smoothed analogs are given in Fig. [7)3. 

FIG. 8. The same as Fig. [?], but now for for a West-O'Donnell power spectrum defined by fc_ = 0.82(o;/c) and k + = 1.97(w/c). 
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FIG. 9. The same as Fig. pk, but now for an angle of incidence of 9q = 5°. 



FIG. 10. The contribution to the mean differential reflection coefficient from the incoherent component of the scattered light 
{dR P /d9) incoh as a function of the scattering angle 9 when a p-polarized plane wave of wavelength A = 633 nm is normally 
incident (Bo = 0°) on the scattering system depicted in Fig. |l|. The dielectric function of the film is e{lo) = 5.6644 + i0.01, and its 
mean thickness is d — 380 nm. The surface profile function £(xi), is characterized by a Gaussian surface height distribution with 
an RMS- height, a = 30 nm and a Gaussian surface-height autocorrelation function with a 1/e-correlation length of a = 100 nm. 
The scattering system supports six satellite peaks (see Ref. |n]]), but two of them are in the non-radiative part of the spectrum. 
The real peaks should appear at 07j 2 , = ±17.5° and 9% 3 . = ±46.1°. Notice that only the smoothed data are given. 



FIG. 11. The same as Fig. 10, but now for a West-O'Donnell power spectrum defined by fc_ = 1.61(lu/c) and k+ = 2.76(cj/c) 
for an angle of incidence 8o = 0° (a) and 9o = 5° (b). The real satellite peaks are predicted to appear at 9 ± = ±17.5° (9q = 0°) 
and 8 ± = 12.3°, —22.8° (9o = 5°). Only the smoothed data are presented. 
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Figure 8 
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